remove(list = ls())
library("lavaan")
library("lavaan.survey")
library("dplyr")
library("plyr")
library("foreign")
library("semTools")
library("xtable")
library("survey")
library("tableone")
library("readstata13")
library("tidyr")

completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

location <- "/Users/dillonlaaker/Box/Stability/"
setwd(paste(location, "data/datasets", sep = ""))
liss <- read.dta13("liss_data.dta")
bes <- read.dta13("bes_data.dta")

####################
####Tables A6-A9####
####################

##########
###LISS###
##########

liss$missing[liss$missing8==1] <- NA #Removes respondents who were not in wave 1
liss <- completeFun(liss, "missing")

liss <- liss %>% mutate(complete = ifelse(missing==0, 1, 0)) #Creates binary variable equal to 1 if respondent completed all waves 0 otherwise

liss <- liss %>% mutate(complete1 = ifelse(missing==7, 0, 1)) #Creates binary variable equal to 1 if respondent completed all waves 0 otherwise

#liss <- liss %>% mutate(complete1 = ifelse(missing==0, 1, ifelse(missing==8, 0, NA)))#Creates binary variable equal to 1 if respondent complete all waves and 0 if the respondent all completed 1 wave.

liss$complete <- ifelse(liss$missing==0, 1, 0) #Creates binary variable equal to 1 if respondent completed all waves 0 otherwise
liss$complete1 <- ifelse(liss$missing==0, 1, ifelse(liss$missing==8, 0, NA)) #Creates binary variable equal to 1 if respondent complete all waves and 0 if the respondent all completed 1 wave.

liss <- liss %>% select(immigration_17, immigration_16, immigration_14, immigration_13, immigration_12, immigration_11, immigration_10, immigration_9, immigration_8, imm1_8, imm1_9, imm1_10, imm1_11, imm1_12, imm1_13, imm1_14, imm1_16, imm1_17, imm2_8, imm2_9, imm2_10, imm2_11, imm2_12, imm2_13, imm2_14, imm2_16, imm2_17, imm3_8, imm3_9, imm3_10, imm3_11, imm3_12, imm3_13, imm3_14, imm3_16, imm3_17, imm4_8, imm4_9, imm4_10, imm4_11, imm4_12, imm4_13, imm4_14, imm4_16, imm4_17, imm5_8, imm5_9, imm5_10, imm5_11, imm5_12, imm5_13, imm5_14, imm5_16, imm5_17, imm6_8, imm6_9, imm6_10, imm6_11, imm6_12, imm6_13, imm6_14, imm6_16, imm6_17, newsinterest8, polinterest8, ideology8, econsatisfied8, econconfidence8, education8, age8, female8, complete, complete1)

xvars1 <- c("immigration_8", "newsinterest8", "age8", "econconfidence8", "econsatisfied8", "polinterest8", "ideology8", "female8", "education8")
labels <- c("Immigration", "News Interest", "Age", "Econ. Confi.", "Econ. Satis.", "Pol. Interest", "Ideology", "Female", "Education") 
liss <- completeFun(liss, xvars1)

##Comparing between respondents who completed all waves and those who did not
nwdata <- svydesign(ids=~1, data=liss)
table <- svyCreateTableOne(vars=xvars1, strata="complete", data=nwdata, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA7a.tex", sep = ""))

model <- glm(complete ~ age8 + I(age8^2) + newsinterest8 + econconfidence8 + econsatisfied8 + polinterest8 + ideology8 + female8 + education8, data=liss, family=binomial(link = "probit"))
ps <- predict(model, type="response")
liss <- liss %>% mutate(weights = ifelse(complete==1, 1/ps, 1/(1-ps)))

wdata1 <- svydesign(ids=~1, data=liss, weights=liss$weights)
table <- svyCreateTableOne(vars=xvars1, strata="complete", data=wdata1, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA7c.tex", sep = ""))
lissa <- liss

##Comparing between respondents who completed all waves and those who completed one wave
liss <- completeFun(liss, "complete1")
nwdata <- svydesign(ids=~1, data=liss)
table <- svyCreateTableOne(vars=xvars1, strata="complete1", data=nwdata, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA9a.tex", sep = ""))

model <- glm(complete1 ~ age8 + I(age8^2) + newsinterest8 + econconfidence8 + econsatisfied8 + polinterest8 + ideology8 + female8 + education8, data=liss, family=binomial(link = "probit"))
ps <- predict(model, type="response")
liss <- liss %>% mutate(weights = ifelse(complete1==1, 1/ps, 1/(1-ps)))

wdata1 <- svydesign(ids=~1, data=liss, weights=liss$weights)
table <- svyCreateTableOne(vars=xvars1, strata="complete1", data=wdata1, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA9c.tex", sep = ""))
lissb <- liss

attach(lissa)
lissa <- data.frame(immigration_17, immigration_16, immigration_14, immigration_13, immigration_12, immigration_11, immigration_10, immigration_9, immigration_8, age8, weights, imm1_8, imm1_9, imm1_10, imm1_11, imm1_12, imm1_13, imm1_14, imm1_16, imm1_17, imm2_8, imm2_9, imm2_10, imm2_11, imm2_12, imm2_13, imm2_14, imm2_16, imm2_17, imm3_8, imm3_9, imm3_10, imm3_11, imm3_12, imm3_13, imm3_14, imm3_16, imm3_17, imm4_8, imm4_9, imm4_10, imm4_11, imm4_12, imm4_13, imm4_14, imm4_16, imm4_17, imm5_8, imm5_9, imm5_10, imm5_11, imm5_12, imm5_13, imm5_14, imm5_16, imm5_17, imm6_8, imm6_9, imm6_10, imm6_11, imm6_12, imm6_13, imm6_14, imm6_16, imm6_17)
detach(lissa)

attach(lissb)
lissb <- data.frame(immigration_17, immigration_16, immigration_14, immigration_13, immigration_12, immigration_11, immigration_10, immigration_9, immigration_8, age8, weights, imm1_8, imm1_9, imm1_10, imm1_11, imm1_12, imm1_13, imm1_14, imm1_16, imm1_17, imm2_8, imm2_9, imm2_10, imm2_11, imm2_12, imm2_13, imm2_14, imm2_16, imm2_17, imm3_8, imm3_9, imm3_10, imm3_11, imm3_12, imm3_13, imm3_14, imm3_16, imm3_17, imm4_8, imm4_9, imm4_10, imm4_11, imm4_12, imm4_13, imm4_14, imm4_16, imm4_17, imm5_8, imm5_9, imm5_10, imm5_11, imm5_12, imm5_13, imm5_14, imm5_16, imm5_17, imm6_8, imm6_9, imm6_10, imm6_11, imm6_12, imm6_13, imm6_14, imm6_16, imm6_17)
detach(lissb)
lissa <- na.omit(lissa)
lissb <- na.omit(lissb)


#########
###BES###
#########

bes$missing[bes$missing1==1] <- NA
bes <- completeFun(bes, "missing")

bes <- bes %>% mutate(complete = ifelse(missing==0, 1, 0)) #Creates binary variable equal to 1 if respondent completed all waves 0 otherwise
bes <- bes %>% mutate(complete1 = ifelse(missing==7, 0, 1))

#bes <- bes %>% mutate(complete1 = ifelse(missing==0, 1, ifelse(missing==7, 0, NA)))#Creates binary variable equal to 1 if respondent complete all waves and 0 if the respondent only completed 1 wave.

bes <- bes %>% select(immecon1, immecon2, immecon3, immecon4, immecon7, immecon8, immecon10, immecon11, immcul1, immcul2, immcul3, immcul4, immcul7, immcul8, immcul10, immcul11, immwelfare1, immwelfare2, immwelfare3, immwelfare4, immwelfare7, immwelfare8, immwelfare10, immwelfare11, imm1, imm2, imm3, imm4, imm7, imm8, imm10, imm11, ideology1, education1, polattention1, econpersretro1, econgenretro1, female, ethnicity, income_house, complete, complete1, age)

bes <- completeFun(bes, c("imm1", "polattention1", "age", "econpersretro1", "econgenretro1", "ideology1", "female", "income_house", "ethnicity", "education1"))
xvars1 <- c("imm1", "polattention1", "age", "econpersretro1", "econgenretro1", "ideology1", "female", "income_house", "ethnicity", "education1")
labels <- c("Immigration", "Political Attention", "Age", "Self Econ", "General Econ", "Ideology", "Female", "Income", "Ethnicity", "Education")


##Comparing between respondents who completed all waves and those who completed one wave
nwdata <- svydesign(ids=~1, data=bes)
table <- svyCreateTableOne(vars=xvars1, strata="complete", data=nwdata, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA7b.tex", sep = ""))

model <- glm(complete ~ age + I(age^2) + polattention1 + econgenretro1 + econpersretro1 + ideology1 + female + income_house + ethnicity + education1, data=bes, family=binomial(link = "probit"))
ps<- predict(model, type="response")
bes <- bes %>% mutate(weights = ifelse(complete==1, 1/(ps), 1/(1-ps)))

wdata1 <- svydesign(ids=~1, data=bes, weights=bes$weights)
table <- svyCreateTableOne(vars=xvars1, strata="complete", data=wdata1, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA7d.tex", sep = ""))
besa <- bes

##Comparing between respondents who completed all waves and those who completed one wave
bes <- completeFun(bes, "complete1")
nwdata <- svydesign(ids=~1, data=bes)
table <- svyCreateTableOne(vars=xvars1, strata="complete1", data=nwdata, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA9b.tex", sep = ""))

model <- glm(complete1 ~ age + I(age^2) + polattention1 + econgenretro1 + econpersretro1 + ideology1 + female + income_house + ethnicity + education1, data=bes, family=binomial(link = "probit"))
ps<- predict(model, type="response")
bes <- bes %>% mutate(weights = ifelse(complete1==1, 1/(ps), 1/(1-ps)))

wdata1 <- svydesign(ids=~1, data=bes, weights=bes$weights)
table <- svyCreateTableOne(vars=xvars1, strata="complete1", data=wdata1, test=FALSE)
table <- print(table, smd=TRUE)
table <- print(table[-1,], smd=TRUE)
rownames(table) <- labels
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=TRUE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "Draft/tables/tableA9d.tex", sep = ""))
besb <- bes


attach(besa)
besa <- data.frame(immecon1, immecon2, immecon3, immecon4, immecon7, immecon8, immecon10, immecon11, immcul1, immcul2, immcul3, immcul4, immcul7, immcul8, immcul10, immcul11, immwelfare1, immwelfare2, immwelfare3, immwelfare4, immwelfare7, immwelfare8, immwelfare10, immwelfare11, imm1, imm2, imm3, imm4, imm7, imm8, imm10, imm11, weights, age)
detach(besa)
attach(besb)
besb <- data.frame(immecon1, immecon2, immecon3, immecon4, immecon7, immecon8, immecon10, immecon11, immcul1, immcul2, immcul3, immcul4, immcul7, immcul8, immcul10, immcul11, immwelfare1, immwelfare2, immwelfare3, immwelfare4, immwelfare7, immwelfare8, immwelfare10, immwelfare11, imm1, imm2, imm3, imm4, imm7, imm8, imm10, imm11, weights, age)
detach(besb)
besa <- na.omit(besa)
besb <- na.omit(besb)



#########################
########Table A10########
#########################
#Uses specification that compares respondents who completed all waves and those who did not

model1 <- "imm8 =~ imm5_8 + imm2_8 + imm3_8 + imm4_8 + imm1_8 + imm6_8
imm9 =~ imm5_9 + imm2_9 + imm3_9 + imm4_9 + imm1_9 + imm6_9
imm10 =~ imm5_10 + imm2_10 + imm3_10 + imm4_10 + imm1_10 + imm6_10
imm11 =~ imm5_11 + imm2_11 + imm3_11 + imm4_11 + imm1_11 + imm6_11
imm12 =~ imm5_12 + imm2_12 + imm3_12 + imm4_12 + imm1_12 + imm6_12
imm13 =~ imm5_13 + imm2_13 + imm3_13 + imm4_13 + imm1_13 + imm6_13
imm14 =~ imm5_14 + imm2_14 + imm3_14 + imm4_14 + imm1_14 + imm6_14
imm16 =~ imm5_16 + imm2_16 + imm3_16 + imm4_16 + imm1_16 + imm6_16
imm17 =~ imm5_17 + imm2_17 + imm3_17 + imm4_17 + imm1_17 + imm6_17
imm9 ~ imm8
imm10 ~ imm9
imm11 ~ imm10
imm12 ~ imm11
imm13 ~ imm12
imm14 ~ imm13
imm16 ~ imm14
imm17 ~ imm16"
fit1 <- sem(model1, se="robust.huber.white", test="Satorra-Bentler", data=lissa)
wliss <- svydesign(ids = ~1, data = lissa, weights=lissa$weights)
fit1 <- lavaan.survey(fit1, wliss)
reliabl1 <- reliability(fit1)
est <- fitMeasures(fit1, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model1results <- parameterEstimates(fit1)
model1results <- model1results[c(55,56,57,58,59,60,61,62),c(4,5)]
model1results <- model1results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 8)
par2 <- rep(")", 8)
model1results$se <- paste(par1, model1results$se, par2, sep="")
model1results$est <- paste(model1results$est, model1results$se)
model1results$se <- NULL
coefnames <- c("$\\beta_{1,2}$", "$\\beta_{2,3}$", "$\\beta_{3,4}$", "$\\beta_{4,5}$", "$\\beta_{5,6}$", "$\\beta_{6,7}$", "$\\beta_{7,8}$", "$\\beta_{8,9}$")
model1results <- data.frame(coefnames, model1results)
coefnames1 <- c("$\\chi^{2}$", "df", "p-value", "CFI", "TLI", "RMSEA", "SRMSR", "AIC")
fitmeasures <- data.frame(est)
est <- fitmeasures %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit1)
obs <- data.frame("N", obs)
fitmeasures <- data.frame(coefnames1, est)
colnames(fitmeasures) <- c("coefnames", "est")
colnames(obs) <- c("coefnames", "est")
results1 <- rbind(model1results, fitmeasures, obs)

model2 <- "imm8 =~ imm5_8 + imm2_8 + imm3_8 + imm4_8 + imm1_8 + imm6_8
imm9 =~ imm5_9 + imm2_9 + imm3_9 + imm4_9 + imm1_9 + imm6_9
imm10 =~ imm5_10 + imm2_10 + imm3_10 + imm4_10 + imm1_10 + imm6_10
imm11 =~ imm5_11 + imm2_11 + imm3_11 + imm4_11 + imm1_11 + imm6_11
imm12 =~ imm5_12 + imm2_12 + imm3_12 + imm4_12 + imm1_12 + imm6_12
imm13 =~ imm5_13 + imm2_13 + imm3_13 + imm4_13 + imm1_13 + imm6_13
imm14 =~ imm5_14 + imm2_14 + imm3_14 + imm4_14 + imm1_14 + imm6_14
imm16 =~ imm5_16 + imm2_16 + imm3_16 + imm4_16 + imm1_16 + imm6_16
imm17 =~ imm5_17 + imm2_17 + imm3_17 + imm4_17 + imm1_17 + imm6_17
imm1_8 ~~ imm1_9 + imm1_10 + imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_9 ~~ imm1_10 + imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_10 ~~ imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_11 ~~ imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_12 ~~ imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_13 ~~ imm1_14 + imm1_16 + imm1_17
imm1_14 ~~ imm1_16 + imm1_17
imm1_16 ~~ imm1_17
imm2_8 ~~ imm2_9 + imm2_10 + imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_9 ~~ imm2_10 + imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_10 ~~ imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_11 ~~ imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_12 ~~ imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_13 ~~ imm2_14 + imm2_16 + imm2_17
imm2_14 ~~ imm2_16 + imm2_17
imm2_16 ~~ imm2_17
imm3_8 ~~ imm3_9 + imm3_10 + imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_9 ~~ imm3_10 + imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_10 ~~ imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_11 ~~ imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_12 ~~ imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_13 ~~ imm3_14 + imm3_16 + imm3_17
imm3_14 ~~ imm3_16 + imm3_17
imm3_16 ~~ imm3_17
imm4_8 ~~ imm4_9 + imm4_10 + imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_9 ~~ imm4_10 + imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_10 ~~ imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_11 ~~ imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_12 ~~ imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_13 ~~ imm4_14 + imm4_16 + imm4_17
imm4_14 ~~ imm4_16 + imm4_17
imm4_16 ~~ imm4_17
imm5_8 ~~ imm5_9 + imm5_10 + imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_9 ~~ imm5_10 + imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_10 ~~ imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_11 ~~ imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_12 ~~ imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_13 ~~ imm5_14 + imm5_16 + imm5_17
imm5_14 ~~ imm5_16 + imm5_17
imm5_16 ~~ imm5_17
imm6_8 ~~ imm6_9 + imm6_10 + imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_9 ~~ imm6_10 + imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_10 ~~ imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_11 ~~ imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_12 ~~ imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_13 ~~ imm6_14 + imm6_16 + imm6_17
imm6_14 ~~ imm6_16 + imm6_17
imm6_16 ~~ imm6_17
imm2_8 ~~ imm4_8
imm5_8 ~~ imm6_8
imm4_8 ~~ imm6_8
imm2_8 ~~ imm3_8
imm2_8 ~~ imm6_8
imm1_9 ~~ imm2_9 + imm3_9 + imm4_9 + imm5_9 + imm6_9
imm2_9 ~~ imm3_9 + imm4_9 + imm5_9 + imm6_9
imm3_9 ~~ imm4_9 + imm5_9 + imm6_9
imm4_9 ~~ imm5_9 + imm6_9
imm5_9 ~~ imm6_9
imm1_10 ~~ imm2_10 + imm3_10 + imm4_10 + imm5_10 + imm6_10
imm2_10 ~~ imm3_10 + imm4_10 + imm5_10 + imm6_10
imm3_10 ~~ imm4_10 + imm5_10 + imm6_10
imm4_10 ~~ imm5_10 + imm6_10
imm5_10 ~~ imm6_10
imm1_11 ~~ imm2_11 + imm3_11 + imm4_11 + imm5_11 + imm6_11
imm2_11 ~~ imm3_11 + imm4_11 + imm5_11 + imm6_11
imm3_11 ~~ imm4_11 + imm5_11 + imm6_11
imm4_11 ~~ imm5_11 + imm6_11
imm5_11 ~~ imm6_11
imm1_12 ~~ imm2_12 + imm3_12 + imm4_12 + imm5_12 + imm6_12
imm2_12 ~~ imm3_12 + imm4_12 + imm5_12 + imm6_12
imm3_12 ~~ imm4_12 + imm5_12 + imm6_12
imm4_12 ~~ imm5_12 + imm6_12
imm5_12 ~~ imm6_12
imm1_13 ~~ imm2_13 + imm3_13 + imm4_13 + imm5_13 + imm6_13
imm2_13 ~~ imm3_13 + imm4_13 + imm5_13 + imm6_13
imm3_13 ~~ imm4_13 + imm5_13 + imm6_13
imm4_13 ~~ imm5_13 + imm6_13
imm5_13 ~~ imm6_13
imm1_14 ~~ imm2_14 + imm3_14 + imm4_14 + imm5_14 + imm6_14
imm2_14 ~~ imm3_14 + imm4_14 + imm5_14 + imm6_14
imm3_14 ~~ imm4_14 + imm5_14 + imm6_14
imm4_14 ~~ imm5_14 + imm6_14
imm5_14 ~~ imm6_14
imm1_16 ~~ imm2_16 + imm3_16 + imm4_16 + imm5_16 + imm6_16
imm2_16 ~~ imm3_16 + imm4_16 + imm5_16 + imm6_16
imm3_16 ~~ imm4_16 + imm5_16 + imm6_16
imm4_16 ~~ imm5_16 + imm6_16
imm5_16 ~~ imm6_16
imm5_17 ~~ imm6_17
imm2_17 ~~ imm3_17
imm2_17 ~~ imm6_17
imm2_17 ~~ imm4_17
imm3_17 ~~ imm6_17
imm5_17 ~~ imm4_17
imm9 ~ imm8
imm10 ~ imm9
imm11 ~ imm10
imm12 ~ imm11
imm13 ~ imm12
imm14 ~ imm13
imm16 ~ imm14
imm17 ~ imm16"

fit2<-sem(model2, se="robust.huber.white", test="Satorra-Bentler", data=lissa)
fit2 <- lavaan.survey(fit2, wliss)
reliabl2 <- reliability(fit2)
est <- fitMeasures(fit2, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model2results <- parameterEstimates(fit2)
model2results <- model2results[c(387,388,389,390,391,392,393,394),c(4,5)]
model2results <- model2results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 8)
par2 <- rep(")", 8)
model2results$se <- paste(par1, model2results$se, par2, sep="")
model2results$est <- paste(model2results$est, model2results$se)
model2results$se <- NULL
fitmeasures2 <- data.frame(est)
fitmeasures2 <- fitmeasures2 %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit2)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results2 <- rbind(model2results, fitmeasures2, obs)

lissresults <- cbind(results1, results2)
colnames(lissresults) <- c("", "Model (1)", "Model (2)")

##BES
model1 <- "imm1 =~ immcul1 + immwelfare1 + immecon1
imm2 =~ immcul2 + immwelfare2 + immecon2
imm3 =~ immcul3 + immwelfare3 + immecon3
imm4 =~ immcul4 + immwelfare4 + immecon4
imm7 =~ immcul7 + immwelfare7 + immecon7
imm8 =~ immcul8 + immwelfare8 + immecon8
imm10 =~ immcul10 + immwelfare10 + immecon10
imm11 =~ immcul11 + immwelfare11 + immecon11
imm2 ~ imm1
imm3 ~ imm2
imm4 ~ imm3
imm7 ~ imm4
imm8 ~ imm7
imm10 ~ imm8
imm11 ~ imm10"
fit1 <- sem(model1, se="robust.huber.white", test="Satorra-Bentler", data=besa)
wbes <- svydesign(ids = ~1, data = besa, weights=besa$weights)
fit1 <- lavaan.survey(fit1, wbes)
reliabb1 <- reliability(fit1)
est <- fitMeasures(fit1, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model1results <- parameterEstimates(fit1)
model1results <- model1results[c(25,26,27,28,29,30,31),c(4,5)]
model1results <- model1results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 7)
par2 <- rep(")", 7)
model1results$se <- paste(par1, model1results$se, par2, sep="")
model1results$est <- paste(model1results$est, model1results$se)
model1results$se <- NULL
model1results <- data.frame(model1results)
fix <- data.frame("")
colnames(fix) <- c("est")
model1results <- rbind(model1results, fix)
fitmeasures <- data.frame(est)
est <- fitmeasures %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit1)
obs <- data.frame(obs)
colnames(obs) <- c("est")
fitmeasures <- data.frame(est)
results1 <- rbind(model1results, fitmeasures, obs)

model2 <- "imm1 =~ immcul1 + immwelfare1 + immecon1
imm2 =~ immcul2 + immwelfare2 + immecon2
imm3 =~ immcul3 + immwelfare3 + immecon3
imm4 =~ immcul4 + immwelfare4 + immecon4
imm7 =~ immcul7 + immwelfare7 + immecon7
imm8 =~ immcul8 + immwelfare8 + immecon8
imm10 =~ immcul10 + immwelfare10 + immecon10
imm11 =~ immcul11 + immwelfare11 + immecon11
immecon1 ~~ immecon2 + immecon3 + immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon2 ~~ immecon3 + immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon3 ~~ immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon4 ~~ immecon7 + immecon8 + immecon10 + immecon11
immecon7 ~~ immecon8 + immecon10 + immecon11
immecon8 ~~ immecon10 + immecon11
immecon10 ~~ immecon11
immcul1 ~~ immcul2 + immcul3 + immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul2 ~~ immcul3 + immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul3 ~~ immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul4 ~~ immcul7 + immcul8 + immcul10 + immcul11
immcul7 ~~ immcul8 + immcul10 + immcul11
immcul8 ~~ immcul10 + immcul11
immcul10 ~~ immcul11
immwelfare1 ~~ immwelfare2 + immwelfare3 + immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare2 ~~ immwelfare3 + immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare3 ~~ immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare4 ~~ immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare7 ~~ immwelfare8 + immwelfare10 + immwelfare11
immwelfare8 ~~ immwelfare10 + immwelfare11
immwelfare10 ~~ immwelfare11
immcul1 ~~ immwelfare1
immcul2 ~~ immwelfare2
immcul3 ~~ immwelfare3
immecon4 ~~ immcul4 + immwelfare4
immecon7 ~~ immcul7 + immwelfare7
immecon8 ~~ immcul8 + immwelfare8
immecon10 ~~ immcul10 + immwelfare10
immecon11 ~~ immcul11 + immwelfare11
imm2 ~ imm1
imm3 ~ imm2
imm4 ~ imm3
imm7 ~ imm4
imm8 ~ imm7
imm10 ~ imm8
imm11 ~ imm10"
fit2 <- sem(model2, se="robust.huber.white", test="Satorra-Bentler", data=besa)
fit2 <- lavaan.survey(fit2, wbes)
reliabb2 <- reliability(fit2)
est <- fitMeasures(fit2, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model2results <- parameterEstimates(fit2)
model2results <- model2results[c(122,123,124,125,126,127,128),c(4,5)]
model2results <- model2results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 7)
par2 <- rep(")", 7)
model2results$se <- paste(par1, model2results$se, par2, sep="")
model2results$est <- paste(model2results$est, model2results$se)
model2results$se <- NULL
fix <- data.frame("")
colnames(fix) <- c("est")
model2results <- rbind(model2results, fix)
fitmeasures2 <- data.frame(est)
fitmeasures2 <- fitmeasures2 %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit2)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results2 <- rbind(model2results, fitmeasures2, obs)

besresults <- cbind(results1, results2)
colnames(besresults) <- c("Model (1)", "Model (2)")

fix <- data.frame(c(rep("", times=17)))
table <- cbind(lissresults, fix, besresults)
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
comment          <- list()
comment$pos      <- list()
comment$pos[[1]] <- 1
comment$pos[[2]] <- 2
comment$pos[[3]] <- 3
comment$pos[[4]] <- 4
comment$pos[[5]] <- 5
comment$pos[[6]] <- 6
comment$pos[[7]] <- 7
comment$pos[[8]] <- 8
comment$command  <- c(rep("[1.5ex] \n", times=7), "\\midrule \n")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE, add.to.row = comment,   hline.after = NULL, file=paste(location, "Draft/tables/tableA9.tex", sep = ""))



#########################
########Table A11########
#########################
#Uses specification that compares respondents who completed all waves and those who only completed 1 wave

model1 <- "imm8 =~ imm5_8 + imm2_8 + imm3_8 + imm4_8 + imm1_8 + imm6_8
imm9 =~ imm5_9 + imm2_9 + imm3_9 + imm4_9 + imm1_9 + imm6_9
imm10 =~ imm5_10 + imm2_10 + imm3_10 + imm4_10 + imm1_10 + imm6_10
imm11 =~ imm5_11 + imm2_11 + imm3_11 + imm4_11 + imm1_11 + imm6_11
imm12 =~ imm5_12 + imm2_12 + imm3_12 + imm4_12 + imm1_12 + imm6_12
imm13 =~ imm5_13 + imm2_13 + imm3_13 + imm4_13 + imm1_13 + imm6_13
imm14 =~ imm5_14 + imm2_14 + imm3_14 + imm4_14 + imm1_14 + imm6_14
imm16 =~ imm5_16 + imm2_16 + imm3_16 + imm4_16 + imm1_16 + imm6_16
imm17 =~ imm5_17 + imm2_17 + imm3_17 + imm4_17 + imm1_17 + imm6_17
imm9 ~ imm8
imm10 ~ imm9
imm11 ~ imm10
imm12 ~ imm11
imm13 ~ imm12
imm14 ~ imm13
imm16 ~ imm14
imm17 ~ imm16"
fit1 <- sem(model1, se="robust.huber.white", test="Satorra-Bentler", data=lissb)
wliss <- svydesign(ids = ~1, data = lissb, weights=lissa$weights)
fit1 <- lavaan.survey(fit1, wliss)
reliabl1 <- reliability(fit1)
est <- fitMeasures(fit1, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model1results <- parameterEstimates(fit1)
model1results <- model1results[c(55,56,57,58,59,60,61,62),c(4,5)]
model1results <- model1results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 8)
par2 <- rep(")", 8)
model1results$se <- paste(par1, model1results$se, par2, sep="")
model1results$est <- paste(model1results$est, model1results$se)
model1results$se <- NULL
coefnames <- c("$\\beta_{1,2}$", "$\\beta_{2,3}$", "$\\beta_{3,4}$", "$\\beta_{4,5}$", "$\\beta_{5,6}$", "$\\beta_{6,7}$", "$\\beta_{7,8}$", "$\\beta_{8,9}$")
model1results <- data.frame(coefnames, model1results)
coefnames1 <- c("$\\chi^{2}$", "df", "p-value", "CFI", "TLI", "RMSEA", "SRMSR", "AIC")
fitmeasures <- data.frame(est)
est <- fitmeasures %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit1)
obs <- data.frame("N", obs)
fitmeasures <- data.frame(coefnames1, est)
colnames(fitmeasures) <- c("coefnames", "est")
colnames(obs) <- c("coefnames", "est")
results1 <- rbind(model1results, fitmeasures, obs)

model2 <- "imm8 =~ imm5_8 + imm2_8 + imm3_8 + imm4_8 + imm1_8 + imm6_8
imm9 =~ imm5_9 + imm2_9 + imm3_9 + imm4_9 + imm1_9 + imm6_9
imm10 =~ imm5_10 + imm2_10 + imm3_10 + imm4_10 + imm1_10 + imm6_10
imm11 =~ imm5_11 + imm2_11 + imm3_11 + imm4_11 + imm1_11 + imm6_11
imm12 =~ imm5_12 + imm2_12 + imm3_12 + imm4_12 + imm1_12 + imm6_12
imm13 =~ imm5_13 + imm2_13 + imm3_13 + imm4_13 + imm1_13 + imm6_13
imm14 =~ imm5_14 + imm2_14 + imm3_14 + imm4_14 + imm1_14 + imm6_14
imm16 =~ imm5_16 + imm2_16 + imm3_16 + imm4_16 + imm1_16 + imm6_16
imm17 =~ imm5_17 + imm2_17 + imm3_17 + imm4_17 + imm1_17 + imm6_17
imm1_8 ~~ imm1_9 + imm1_10 + imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_9 ~~ imm1_10 + imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_10 ~~ imm1_11 + imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_11 ~~ imm1_12 + imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_12 ~~ imm1_13 + imm1_14 + imm1_16 + imm1_17
imm1_13 ~~ imm1_14 + imm1_16 + imm1_17
imm1_14 ~~ imm1_16 + imm1_17
imm1_16 ~~ imm1_17
imm2_8 ~~ imm2_9 + imm2_10 + imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_9 ~~ imm2_10 + imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_10 ~~ imm2_11 + imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_11 ~~ imm2_12 + imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_12 ~~ imm2_13 + imm2_14 + imm2_16 + imm2_17
imm2_13 ~~ imm2_14 + imm2_16 + imm2_17
imm2_14 ~~ imm2_16 + imm2_17
imm2_16 ~~ imm2_17
imm3_8 ~~ imm3_9 + imm3_10 + imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_9 ~~ imm3_10 + imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_10 ~~ imm3_11 + imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_11 ~~ imm3_12 + imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_12 ~~ imm3_13 + imm3_14 + imm3_16 + imm3_17
imm3_13 ~~ imm3_14 + imm3_16 + imm3_17
imm3_14 ~~ imm3_16 + imm3_17
imm3_16 ~~ imm3_17
imm4_8 ~~ imm4_9 + imm4_10 + imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_9 ~~ imm4_10 + imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_10 ~~ imm4_11 + imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_11 ~~ imm4_12 + imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_12 ~~ imm4_13 + imm4_14 + imm4_16 + imm4_17
imm4_13 ~~ imm4_14 + imm4_16 + imm4_17
imm4_14 ~~ imm4_16 + imm4_17
imm4_16 ~~ imm4_17
imm5_8 ~~ imm5_9 + imm5_10 + imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_9 ~~ imm5_10 + imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_10 ~~ imm5_11 + imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_11 ~~ imm5_12 + imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_12 ~~ imm5_13 + imm5_14 + imm5_16 + imm5_17
imm5_13 ~~ imm5_14 + imm5_16 + imm5_17
imm5_14 ~~ imm5_16 + imm5_17
imm5_16 ~~ imm5_17
imm6_8 ~~ imm6_9 + imm6_10 + imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_9 ~~ imm6_10 + imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_10 ~~ imm6_11 + imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_11 ~~ imm6_12 + imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_12 ~~ imm6_13 + imm6_14 + imm6_16 + imm6_17
imm6_13 ~~ imm6_14 + imm6_16 + imm6_17
imm6_14 ~~ imm6_16 + imm6_17
imm6_16 ~~ imm6_17
imm2_8 ~~ imm4_8
imm5_8 ~~ imm6_8
imm4_8 ~~ imm6_8
imm2_8 ~~ imm3_8
imm2_8 ~~ imm6_8
imm1_9 ~~ imm2_9 + imm3_9 + imm4_9 + imm5_9 + imm6_9
imm2_9 ~~ imm3_9 + imm4_9 + imm5_9 + imm6_9
imm3_9 ~~ imm4_9 + imm5_9 + imm6_9
imm4_9 ~~ imm5_9 + imm6_9
imm5_9 ~~ imm6_9
imm1_10 ~~ imm2_10 + imm3_10 + imm4_10 + imm5_10 + imm6_10
imm2_10 ~~ imm3_10 + imm4_10 + imm5_10 + imm6_10
imm3_10 ~~ imm4_10 + imm5_10 + imm6_10
imm4_10 ~~ imm5_10 + imm6_10
imm5_10 ~~ imm6_10
imm1_11 ~~ imm2_11 + imm3_11 + imm4_11 + imm5_11 + imm6_11
imm2_11 ~~ imm3_11 + imm4_11 + imm5_11 + imm6_11
imm3_11 ~~ imm4_11 + imm5_11 + imm6_11
imm4_11 ~~ imm5_11 + imm6_11
imm5_11 ~~ imm6_11
imm1_12 ~~ imm2_12 + imm3_12 + imm4_12 + imm5_12 + imm6_12
imm2_12 ~~ imm3_12 + imm4_12 + imm5_12 + imm6_12
imm3_12 ~~ imm4_12 + imm5_12 + imm6_12
imm4_12 ~~ imm5_12 + imm6_12
imm5_12 ~~ imm6_12
imm1_13 ~~ imm2_13 + imm3_13 + imm4_13 + imm5_13 + imm6_13
imm2_13 ~~ imm3_13 + imm4_13 + imm5_13 + imm6_13
imm3_13 ~~ imm4_13 + imm5_13 + imm6_13
imm4_13 ~~ imm5_13 + imm6_13
imm5_13 ~~ imm6_13
imm1_14 ~~ imm2_14 + imm3_14 + imm4_14 + imm5_14 + imm6_14
imm2_14 ~~ imm3_14 + imm4_14 + imm5_14 + imm6_14
imm3_14 ~~ imm4_14 + imm5_14 + imm6_14
imm4_14 ~~ imm5_14 + imm6_14
imm5_14 ~~ imm6_14
imm1_16 ~~ imm2_16 + imm3_16 + imm4_16 + imm5_16 + imm6_16
imm2_16 ~~ imm3_16 + imm4_16 + imm5_16 + imm6_16
imm3_16 ~~ imm4_16 + imm5_16 + imm6_16
imm4_16 ~~ imm5_16 + imm6_16
imm5_16 ~~ imm6_16
imm5_17 ~~ imm6_17
imm2_17 ~~ imm3_17
imm2_17 ~~ imm6_17
imm2_17 ~~ imm4_17
imm3_17 ~~ imm6_17
imm5_17 ~~ imm4_17
imm9 ~ imm8
imm10 ~ imm9
imm11 ~ imm10
imm12 ~ imm11
imm13 ~ imm12
imm14 ~ imm13
imm16 ~ imm14
imm17 ~ imm16"

fit2<-sem(model2, se="robust.huber.white", test="Satorra-Bentler", data=lissb)
fit2 <- lavaan.survey(fit2, wliss)
reliabl2 <- reliability(fit2)
est <- fitMeasures(fit2, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model2results <- parameterEstimates(fit2)
model2results <- model2results[c(387,388,389,390,391,392,393,394),c(4,5)]
model2results <- model2results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 8)
par2 <- rep(")", 8)
model2results$se <- paste(par1, model2results$se, par2, sep="")
model2results$est <- paste(model2results$est, model2results$se)
model2results$se <- NULL
fitmeasures2 <- data.frame(est)
fitmeasures2 <- fitmeasures2 %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit2)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results2 <- rbind(model2results, fitmeasures2, obs)

lissresults <- cbind(results1, results2)
colnames(lissresults) <- c("", "Model (1)", "Model (2)")

##BES
model1 <- "imm1 =~ immcul1 + immwelfare1 + immecon1
imm2 =~ immcul2 + immwelfare2 + immecon2
imm3 =~ immcul3 + immwelfare3 + immecon3
imm4 =~ immcul4 + immwelfare4 + immecon4
imm7 =~ immcul7 + immwelfare7 + immecon7
imm8 =~ immcul8 + immwelfare8 + immecon8
imm10 =~ immcul10 + immwelfare10 + immecon10
imm11 =~ immcul11 + immwelfare11 + immecon11
imm2 ~ imm1
imm3 ~ imm2
imm4 ~ imm3
imm7 ~ imm4
imm8 ~ imm7
imm10 ~ imm8
imm11 ~ imm10"
fit1 <- sem(model1, se="robust.huber.white", test="Satorra-Bentler", data=besb)
wbes <- svydesign(ids = ~1, data = besb, weights=besa$weights)
fit1 <- lavaan.survey(fit1, wbes)
reliabb1 <- reliability(fit1)
est <- fitMeasures(fit1, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model1results <- parameterEstimates(fit1)
model1results <- model1results[c(25,26,27,28,29,30,31),c(4,5)]
model1results <- model1results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 7)
par2 <- rep(")", 7)
model1results$se <- paste(par1, model1results$se, par2, sep="")
model1results$est <- paste(model1results$est, model1results$se)
model1results$se <- NULL
model1results <- data.frame(model1results)
fix <- data.frame("")
colnames(fix) <- c("est")
model1results <- rbind(model1results, fix)
fitmeasures <- data.frame(est)
est <- fitmeasures %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit1)
obs <- data.frame(obs)
colnames(obs) <- c("est")
fitmeasures <- data.frame(est)
results1 <- rbind(model1results, fitmeasures, obs)

model2 <- "imm1 =~ immcul1 + immwelfare1 + immecon1
imm2 =~ immcul2 + immwelfare2 + immecon2
imm3 =~ immcul3 + immwelfare3 + immecon3
imm4 =~ immcul4 + immwelfare4 + immecon4
imm7 =~ immcul7 + immwelfare7 + immecon7
imm8 =~ immcul8 + immwelfare8 + immecon8
imm10 =~ immcul10 + immwelfare10 + immecon10
imm11 =~ immcul11 + immwelfare11 + immecon11
immecon1 ~~ immecon2 + immecon3 + immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon2 ~~ immecon3 + immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon3 ~~ immecon4 + immecon7 + immecon8 + immecon10 + immecon11
immecon4 ~~ immecon7 + immecon8 + immecon10 + immecon11
immecon7 ~~ immecon8 + immecon10 + immecon11
immecon8 ~~ immecon10 + immecon11
immecon10 ~~ immecon11
immcul1 ~~ immcul2 + immcul3 + immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul2 ~~ immcul3 + immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul3 ~~ immcul4 + immcul7 + immcul8 + immcul10 + immcul11
immcul4 ~~ immcul7 + immcul8 + immcul10 + immcul11
immcul7 ~~ immcul8 + immcul10 + immcul11
immcul8 ~~ immcul10 + immcul11
immcul10 ~~ immcul11
immwelfare1 ~~ immwelfare2 + immwelfare3 + immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare2 ~~ immwelfare3 + immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare3 ~~ immwelfare4 + immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare4 ~~ immwelfare7 + immwelfare8 + immwelfare10 + immwelfare11
immwelfare7 ~~ immwelfare8 + immwelfare10 + immwelfare11
immwelfare8 ~~ immwelfare10 + immwelfare11
immwelfare10 ~~ immwelfare11
immcul1 ~~ immwelfare1
immcul2 ~~ immwelfare2
immcul3 ~~ immwelfare3
immecon4 ~~ immcul4 + immwelfare4
immecon7 ~~ immcul7 + immwelfare7
immecon8 ~~ immcul8 + immwelfare8
immecon10 ~~ immcul10 + immwelfare10
immecon11 ~~ immcul11 + immwelfare11
imm2 ~ imm1
imm3 ~ imm2
imm4 ~ imm3
imm7 ~ imm4
imm8 ~ imm7
imm10 ~ imm8
imm11 ~ imm10"
fit2 <- sem(model2, se="robust.huber.white", test="Satorra-Bentler", data=besb)
fit2 <- lavaan.survey(fit2, wbes)
reliabb2 <- reliability(fit2)
est <- fitMeasures(fit2, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
model2results <- parameterEstimates(fit2)
model2results <- model2results[c(122,123,124,125,126,127,128),c(4,5)]
model2results <- model2results %>% mutate(est= format(round(est, 2), nsmall = 2), se = format(round(se, 2), nsmall = 2))
par1 <- rep("(", 7)
par2 <- rep(")", 7)
model2results$se <- paste(par1, model2results$se, par2, sep="")
model2results$est <- paste(model2results$est, model2results$se)
model2results$se <- NULL
fix <- data.frame("")
colnames(fix) <- c("est")
model2results <- rbind(model2results, fix)
fitmeasures2 <- data.frame(est)
fitmeasures2 <- fitmeasures2 %>% mutate(est= format(round(est, 2), nsmall = 2))
obs <- nobs(fit2)
obs <- data.frame(obs)
colnames(obs) <- c("est")
results2 <- rbind(model2results, fitmeasures2, obs)

besresults <- cbind(results1, results2)
colnames(besresults) <- c("Model (1)", "Model (2)")

fix <- data.frame(c(rep("", times=17)))
table <- cbind(lissresults, fix, besresults)
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
comment          <- list()
comment$pos      <- list()
comment$pos[[1]] <- 1
comment$pos[[2]] <- 2
comment$pos[[3]] <- 3
comment$pos[[4]] <- 4
comment$pos[[5]] <- 5
comment$pos[[6]] <- 6
comment$pos[[7]] <- 7
comment$pos[[8]] <- 8
comment$command  <- c(rep("[1.5ex] \n", times=7), "\\midrule \n")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE, add.to.row = comment,   hline.after = NULL, file=paste(location, "Draft/tables/tableA10.tex", sep = ""))

